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ABSTRACT 

We present numerical magnetohydrostatic solutions describing the gravita- 
tionally stratified, bulk equilibrium of cool, dense prominence plasma embedded 
in a near-potential coronal field. These solutions are calculated using the FI- 
NESSE magnetohydrodynamics equilibrium solver and describe the morphologies 
of magnetic field distributions in and around prominences and the cool promi- 
nence plasma that these fields support. The equilibrium condition for this class 
of problem is usually different in distinct subdomains, separated by free bound- 
aries, across which solutions are matched by suitable continuity or jump con- 
ditions describing force balance. We employ our precise finite element elliptic 
solver to calculate solutions not accessible by previous analytical techniques with 
temperature or entropy prescribed as free functions of the magnetic flux func- 
tion, including a range of values of the polytropic index, temperature variations 
mainly across magnetic field lines and photospheric field profiles sheared close to 
the polarity inversion line. Out of the many examples computed here, perhaps 
the most noteworthy is one which reproduces precisely the three-part structure 
often encountered in observations: a cool dense prominence within a cavity/flux 
rope embedded in a hot corona. The stability properties of these new equilibria, 
which may be relevant to solar eruptions, can be determined in the form of a full 
resistive MHD spectrum using a companion hyperbolic stability solver. 



Subject headings: magnetohydrodynamics: Sun, solar prominences, solar mag- 
netic fields, solar corona 
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1. Introduction 

Solar prominences are cool, dense concentrations of plasma suspended by magnetic 
fields about 10,000 km clear of the base of the corona (Tandberg-Hanssen 1995). The promi- 
nence is referred to as a filament when it is observed in absorption against the solar disk. 
Observing, understanding and modeling prominences and the equilibrium and stability prop- 
erties of their magnetic fields is a major topic of solar research. Eventually many quiescent 
prominences erupt abruptly and are expelled through the corona into interplanetary space 
embedded within a coronal mass ejection (CME). CMEs are the most powerful drivers of the 
Sun-Earth connection. Observations (e.g. Munro et al. 1979, Webb & Hundhausen 1987) 
and recent developments in solar magnetohydrodynamics have shown that more than 70% 
of CMEs originate from eruptions of prominences caused by failures of the confinement of 
highly twisted magnetic fields (Low 1996, 2001, Amari et al. 2003a, 2003b; Fan & Gibson 
2003, 2004; Fan & Low 2003; Flyer et al. 2004, 2005, Zhang & Low 2004, 2005). Recently, 
using magneto-friction models, Mackay & van Ballegooijen (2006) considered the formation 
and ejection of flux ropes as the result of the build-up and concentration of axial flux by a 
simple reconnection process. 

The magnetic field is central to the formation and persistence of quiescent prominences 
(Tandberg-Hanssen 1995). It provides support against gravity for the electrically highly 
conducting prominence plasma and also acts as a thermal shield for the cool prominence 
against the million-degree coronal environment. Observations have shown that long quiescent 
prominences form from the remnants of decaying active regions (Gaizauskas et al. 2001), and 
tend to lie in a cavity at the base of a well-formed coronal helmet streamer, running along 
the polarity inversion line of an extensive bipolar region on the photosphere (Tandberg- 
Hanssen 1995). These prominences also tend to be suspended higher in the corona (> 
30, 000 km above the photosphere) than those in active regions, which tend to be shorter in 
length and more contorted in shape by the more complex active region fields. We adopt the 
hypothesis of Low & Hundhausen (1995) that quiescent and active-region prominences share 
the same basic magnetic structure and mechanics although they differ in magnetic intensity, 
structural complexity, structural length scale and evolutionary time scale. 

The long life of a quiescent prominence as a macroscopic structure, from days to a week 
or more, suggests that static equilibrium is a reasonable first approximation to describe the 
prominence (Anzer 1989). We confine our equilibrium study to 2D configurations for two 
reasons. Firstly, prominences are significantly larger in one horizontal direction than the 
other, being about 200,000 km long and about 5,000 km wide, so that their bulk properties 
can be modeled in 2D. They are invariably found parallel to and above photospheric polarity 
inversion lines. Our 2D models describe flux ropes levitating in the solar atmosphere above 
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the photosphere, which is taken as an infinite plane. The topology of the ends of the flux 
rope, which in reality curve down to meet the photosphere in a full 3D configuration, is left 
out of the 2D description. This description captures features of flux ropes the main part of 
whose lengths are suspended in the atmosphere above the photosphere. 

The second reason is that in 3D the MHS force-balance equation becomes highly non- 
linear, and a follow-up stability analysis becomes even more challenging as the effect of 
line-tying at the photospheric end points must be taken into account as well. Furthermore, 
the 3D MHS force-balance equation changes type from elliptic to mixed type, causing the 
character of the boundary value problem to change. The 2D MHS (Grad-Shafranov) equation 
has a linear, elliptic differential operator so that its characteristic curves are all imaginary. 
Thus, like fields <p described by, e.g., Laplace's equation V 2 = or the linear wave equa- 
tion V 2 + k 2 (p = 0, solutions ip of the Grad-Shafranov equation are continuous functions 
uniquely determined by specifying their value, normal derivative or linear combination of 
the two on the boundary. The situation with the 3D MHS problem is very different (Parker 
1979, 1994): the inclusion of three-dimensional variations causes the force-balance equation 
to change type from a simple elliptic equation with purely imaginary characteristics to an 
equation of mixed type with both imaginary and real characteristics. Across the real charac- 
teristic curves field derivatives may not be defined. Since these real characteristics coincide 
with field trajectories, tangential discontinuities are permitted there provided that the total 
pressure |B| 2 /8 7r + p is balanced across the discontinuity. In general, equilibrium solvers 
break down for mixed-type problems and it can been questioned whether the equilibrium 
and perturbation can be meaningfully distinguished in these mixed-type regimes, as has been 
questioned for hyperbolic regimes in mixed-type problems (Goedbloed 2003, 2004). 

The stability properties of these new equilibria, which may be relevant to solar eruptions, 
can be determined in the form of a full resistive MHD spectrum by solving the linearized 
MHD equations using the companion hyperbolic stability solver PHOENIX (Blokland et 
al. 2007a). MHD waves and instabilities control the dynamics of plasma and occur as the 
natural response to global excitation. Measurement of the spectrum of MHD waves gives 
direct information on the internal state of the plasma. This is called MHD spectroscopy 
in analogy with quantum mechanical spectroscopy, which also involves eigenvalue problems 
of linear operators. MHD spectroscopy entails a separate study of the nonlinear static 
equilibrium configuration on the one hand and the various linear wave structures that can 
occur on the other. The first study is the subject of this paper and the spectroscopy will be 
treated in a sequel. 

The paper is organized as follows. We discuss the morphology of prominences and their 
associated magnetic fields in Section [21 The basic magnetohydrostatic equilibrium problem 



-4- 



is introduced in Section [3] before the new models are presented in Section HI We conclude 
with a discussion in Section [51 

2. Prominence magnetic field morphology 

After a lapse of 20 years or so, there has been a renewed interest in the measurement of 
prominence (and coronal) magnetic fields (Judge 1998, Lin et al. 1998, 2000, Lopez Ariste 
& Casini 2002, 2003, Trujillo Bueno 2003). New techniques for observing magnetic fields of 
prominences on the solar disk (Lin et al. 1998, Trujillo Bueno 2003) allow many new pos- 
sibilities for diagnosing prominences, including low-lying ones inaccessible to older methods 
which relied on the prominences being clearly visible above the limb. 

Spectro-polarimetry has shown that a prominence is threaded by a nearly horizontal 
magnetic field whose principal component is along the long axis of the prominence. Promi- 
nences are classified by whether their magnetic fields thread in the same (called normal 
topology) or opposite (inverse topology) compared to the underlying photospheric bipolar 
field. These two classes imply distinct topologies for the coronal magnetic fields around 
prominences. These two basic topologies, normal and inverse, were first captured in the 
models of Kippenhahn & Schluter (1957) and Kuperus & Raadu 1974), respectively. The 
observed principal component of prominence magnetic fields along the long axes of promi- 
nences implies that fields around prominences are characterized by highly sheared or twisted 
magnetic topologies. 

Statistical studies of prominence topologies suggest that inverse prominences are more 
common than normal prominences by a factor of about three, while inverse prominences 
sit higher in the atmosphere than normal prominences (Leroy et al. 1984, 1989, Bommier et 
al. 1994). Most prominences associated with helmet streamers have the inverse configuration 
(Leroy et al. 1984). In the 2D inverse topology, the prominence must sit in a two-flux 
magnetic system, one flux connecting the bipolar magnetic sources in the photosphere below 
and the other forming a rope which embeds the prominence and runs above and parallel 
to the photospheric polarity inversion line. In 3D, single-flux models can produce inverse- 
topologies and helmet streamer structures (e.g. Aulanier & Demoulin 1998, Mackay & van 
Ballegooijen 2005, 2006). The prominence flux rope is seen as a cavity in coronal white light 
observations and as a filament channel in the chromosphere. For the inverse topology with 
the prominence magnetic field in the opposite direction to that of the bipolar photospheric 
field region below, the magnetic flux rope containing the prominence supports a part of the 
prominence weight by current attraction from above (Low & Hundhausen 1995). Low & 
Zhang (2002) suggested that the two loosely defined classes of impulsively and gradually 
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accelerated CMEs, fast and slow CMEs, may be explained by the different interactions 
between expulsion dynamics and magnetic reconnections taking place in the different inverse 
and normal field topologies. 

A solar prominence cavity or flux rope may be idealized as a 2D structure running 
above and parallel to the photospheric polarity inversion line in the ignored direction, the 
x-direction. In the y-z plane, it is represented by a finite region of closed contours of 
the magnetic flux function with cool plasma localized by gravitational stratification at the 
bottom of the rope. The characteristic hydrostatic scale length of the cool prominence 
plasma is a few hundred km while a typical flux rope size is estimated to be 10 Mm or more. 
While details of magnetic field structure in and around prominences have not been precisely 
determined by observations, Low & Hundhausen (1995) give a plausible sketch of how this 
field may be structured in their Figure 1. The flux contours are compressed at the bottom 
by the cool, dense prominence plasma. 

Amari & Aly (1989, 1992) modeled prominences by embedding massive MHS line cur- 
rents and current sheets in twisted 2D force-free magnetic fields for the first time. Then Low 
& Hundhausen (1995) gave an insightful model of the two topologies of quiescent promi- 
nences using line currents. Taking the plane z = to be the photosphere, the magnetic field 
above this plane around the prominence mass, which is formed along a line (y, z) = (0, z%) 
parallel to the x-axis, with Z2 > 0, is given by the magnetic flux function 

^ = ^log[y 2 + (z-zi) 2 ] + ^log 

where c is the speed of light. The 7 p hoto term causes the photospheric field to be bipolar and 
is due to a subphotospheric virtual line current at (y, z) = (0, z\) with z\ < 0. The I pTom 
term represents the line current associated with the prominence itself as well as an image 
current beneath the photosphere which causes this part of the flux function to vanish at the 
photosphere. Thus, the vertical magnetic field B z at the photosphere is entirely unaffected 
by the presence of the filament-image pair, which only changes the horizontal magnetic field 
By. There is no axial magnetic field component B x in this simple model, which is useful 
for investigating the magnetic field topology around the prominence and does not resolve 
the prominence structure itself. For this simple line-current model, the two basic magnetic 
prominence topologies depend on whether J p hoto and J pr om are of the same sign (inverse 
topology) or opposite signs (normal topology). Figure [Q shows an example of each topology. 
We will use this simple model to guide our numerical modeling in later sections. 

Hundhausen & Low (1994) presented an analytical continuation method for calculating 
polytropic MHS equilibrium states in the Cartesian plane. This method was exploited in Low 
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Fig. 1. — The inverse (left) and normal (right) prominence topologies for prominences rep- 
resented by massive line currents. In the inverse topology the line current flows in the same 
direction as the intrinsic photospheric current, while the currents are oppositely directed in 
the normal topology. The Lorentz force is directed upward in both examples, as it must be 
to balance the weight. 

& Hundhausen (1995) 's study of inverse prominences. Whereas Low & Hundhausen (1995) 
represented prominence plasma enhancements by massive line currents and current sheets, 
Low & Zhang (2004) gave analytical solutions describing plasma enhancements of finite 
width by addressing a free-boundary problem. By this approach they were able to describe 
both normal and inverse magnetic topologies. This was made possible by the assumption of 
a circular boundary between the prominence plasma and the ambient corona and by making 
assumptions simplifying the polytropic equation of state. 

In this paper we present classes of prominence-like solutions calculated numerically using 
the FINite Element Solver for Stationary Equilibria (FINESSE, Belien et al. 2002) which 
allows the computation of stationary, gravitationally stratified axisymmetric or cartesian 
configurations. The two codes FINESSE and PHOENIX (Blokland et al. 2007a) have already 
been used in tandem to determine the stability of tokamaks and accretion disks with steady 
flows (Goedbloed et al. 2004a, b, Blokland et al. 2007b) but have not been adapted for the 
study of solar phenomena although the PARIS code (Belien et al. 1997) which treats the 
gravitational stratification of axisymmetric static loops is incorporated in the functionality 
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of FINESSE. Here, we use a suitably improved version of the FINESSE code to extend 
recent analytic work on translationally invariant magnetohydrostatic equilibria with uniform 
gravity in Cartesian geometry, by allowing a free boundary between the prominence flux 
rope and the ambient coronal magnetic field and allowing one to choose the polytropic index 
freely. The various classes of magnetohydrostatic solutions discussed in what follows each 
lead to a second order PDE, which FINESSE solves in weak form using a Picard iteration. 
We implemented the various forms obtained for this PDE under different choices of freely 
chosen flux functions, along with their scalings, as discussed in the Appendix. We restrict our 
numerical calculations to static solutions which are translationally symmetric, all of which 
fall into an elliptic regime where a split between the equilibrium and the perturbations in a 
forthcoming stability analysis can meaningfully be performed. 



3. The magnetohydrostatic problem 

Consider the static-equilibrium model based on the one-fluid ideal hydromagnetic de- 
scription, denoting the magnetic field, plasma pressure and density by B, p and p, respec- 
tively. The balance of forces is described by 

-3- (V x B) x B - Vp — pgz = 0, (2) 

4:71 

assuming a uniform local gravity of acceleration g in the — z Cartesian direction. Then the 
ideal gas law relates the gas pressure p to the gas density p 

V = ^ P T, (3) 
fj, 

where ks is Boltzmann's constant and p is the mean particle mass for a fully ionized 
(monatomic) hydrogen plasma. The solenoidal condition 



V • B = 0, (4) 

closes the set of equations to determine p, p, and B. To keep the physical problem simple, 
we avoid the complication of a full energy equation by applying in turn two assumptions: 
(1) that the plasma temperature is a flux function T = T(ip) (including the isothermal case 
T = T a constant), and (2) the polytropic case where the entropy s = p/ p 1 = s(ip) is 
a flux function. (The polytropic case with 7 equal to the ratio of specific heats describes 
an adiabatic gas while the isothermal case eliminates heat transport.) The value of the 
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polytropic index, denoted here by 7, is not well known in the corona. Parker (1962) gives 
7 w 1.1 — 1.2 as a guide and this seems to hold 40 years after publication. Any value of 7 
less than 5/3, the adiabatic lapse value, is of potential physical interest. The 7 < 1 cases 
extend the work of Low & Hundhausen (1995) and Low & Zhang (2004) in which 7 is of the 
form 7 = n/(n + l) where n is an integer which is positive in their cases. 

In the 2D case with axial symmetry (d/dx = 0), the magnetic field satisfying equa- 
tion P| has the form 



x3 if, \ dip(y,z) dip(y,z)\ 

where f(y, z) is the axial magnetic field component and ip(y, z) is the magnetic flux function, 
whose isosurfaces contain the magnetic field trajectories. By standard theory (Low 1975), 
the momentum equation (T2J) reduces to 



^ m m+*mA = 0, ( a) 



dp(ijj, z) 



dz 



+ P9 = 0. (7) 

t/i=const 



Here f(y,z) = f(ip) must be a strict function of ip from the projection of the momentum 
equation on the symmetry axis. Here we take the partial z derivatives while keeping ijj 
constant. Equation (j7|) is solved by 



p(ip, z) = px(V») - / p(ip, z')gdz' 





Details of the method used to solve Equation (jHJ) numerically are given in the Appendix. 
For simplicity a Cartesian grid with exponentially varying gridpoint spacing is used to cap- 
ture the detailed structure near the symmetry axis. Thus gridpoints are packed around the 
origin while the domain can extend to 50 units or more, at which distances the field is rather 
insensitive to the details of the structure near the origin. 

The solutions are calculated numerically, solving in both potential and non-potential 
regions simultaneously, exploiting the freedom to define Po(ip) and f 2 {ip) arbitrarily. Piece- 
wise definitions of po and f 2 allow one to calculate a full magnetohydrostatic equilibrium 
embedded in a potential field at once. 
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The mapping from the Cartesian space (y, z) onto the (if), z) space is not one-one in 
general because there may be more than one field line of the same value if) = if) c at a given 
height z. The field lines in the inverse configuration typically have two points at a given 
height, or one in the half-domain y > in which the problem is solved in practice. Then 
under symmetry about the z-axis, single functions p(tp, z) and p(if>, z) characterize the gas 
pressure and density fully. Field lines of the normal topology, on the other hand, may have 
four (two in y > 0) distinct points at the same height. Moreover, the same constant value 
of if) may be found at two field lines: one anchored to the base of the atmosphere and one 
within the levitating flux rope. Even with symmetry about the z-axis, the two points at 
the same height are generally complicated by having distinct gas pressures while they are 
mapped onto the same point (if), z). In such a situation, p(ifi, z) and p(if>, z) are multivalued 
functions of if) (Low & Zhang 2004). 

In general, the Grad-Shafranov equation can assume as many different forms as there 
are distinct regions of physical space, each occupied by a different branch of the functions 
p(if>, z) and p(if), z). The distinct regions are separated by free boundaries, the loci of which 
are unknown and must be calculated simultaneously with the solution ip. The free boundaries 
are governed by continuity or jump conditions imposed by force balance at these boundaries. 
Free boundaries also separate regions where the functional forms of p(if>, z) and p(if>, z) differ, 
e.g., a region of potential or force-free magnetic field separated from a region of full MHS 
force balance. 



Before describing models with MHS force balance we discuss force-free flux-rope solu- 
tions. The force- free case is the one with Po(ip) = constant so that all Maxwell stresses are 
contained within the magnetic field and all electric currents are aligned with the magnetic 
field. The simplest case, the potential case f 2 (if>) = cannot have closed magnetic flux 
contours by the maximum principle for the Laplacian differential operator. Electric currents 
are necessary for a flux rope to be present. The current associated with the magnetic field 
of Equation ([5]) is described by 



4. 



The models 



4.1. 



Nonlinear force-free fields 




(9) 



and the linear case f{if>) = cx-Qip, for some constant ojq, is seen to yield the usual Helmholtz 
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equation for linear force- free fields (V 2 + aQ)ip = - see also Equation ([6]). The general condi- 
tion for the nonlinear field to be force-free is that Equation ([6]) is satisfied with dp(y, z)/dz = 
0, in which case j = a(y, ^)B with the nonlinear force-free parameter a(y,z) = df(ip)/dip. 
The field-aligned currents are clearly associated with axial magnetic flux B x which is respon- 
sible for twist and shear in the magnetic field in these 2D solutions. 

If we insist that the ambient magnetic field be near-potential then the current distribu- 
tion must be concentrated in the vicinity of the flux rope. Therefore, in our examples f 2 (ip) 
is non-zero only in regions bounded by a chosen closed flux contour (here ip = 0.5). The 
observed increase of magnetic field strength and emission intensity to maxima at the top of 
a prominence (Rust 1967, House & Smart 1982, Leroy et al. 1989) suggest forms of f(ip) 
and Po{tp) with maxima at the flux rope center. In Low & Hundhausen (1995) and Low & 
Zhang (2004), f 2 is of the form 1/2 — ip, taking ip — at the magnetic axis at the center 
of the flux rope and ip strictly positive elsewhere. This profile is shown in Figure [2] as curve 
a. Outside the contour ip = 0.5, f 2 (ip) = 0. Figure [2] also shows several profiles for f 2 (ip). 
Curves b, c and d are all of the form 

f(iP) = \(l-2iP) 2 (l + niP), (10) 

with n = 2, 3 and —2 respectively. There functions are the Hermite polynomials with double 
roots at ip = 0.5 and satisfying the conditions at ip = f 2 = 0.5 and df 2 /dip = —1, —1/2 and 
—3 respectively. The double root at ip — 0.5 allows a smooth transition from the current- 
carrying flux rope to the current-free ambient field. The corresponding flux contours for the 
four cases are shown in Figure [3] and the angle of the magnetic field with respect to the 
magnetic axis is shown for each case as functions of space in Figure HI The angle must be 
zero at ip = since B y and B z must vanish there. Outside the ip = 0.5 contour, the angle 
must be 90° since B x = there. Therefore in each case the angle must vary from 0° to 90° 
within the flux rope. The value of angle that an observer would find from one of our models 
depends on the nature of this variation from 0° to 90° and on the distribution of plasma 
within the flux rope. Typical measured angles are of order 30° (Leroy et al. 1984, 1989, 
Bommier et al. 1994). 

Magnetic field models corresponding to the four f 2 profiles of Figure [2] are shown in 
Figure [3j The boundary conditions are given by Equation describing the solutions shown 
in Figure [TJ These boundary conditions are motivated by the fact that a (2D) flux rope 
resembles from a distance a line current and because of the relationship of this work to 
Low & Hundhausen (1995), whose models the present work extends. The I pTO m term has 
no contribution to the photospheric boundary condition while its influence on the flux rope 
solutions when applied at the distant lateral and top boundaries is found to be effectively 
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1* 

Fig. 2. — Examples of f 2 {ip) functions applied in this paper. The quantity f(tp) is the axial 
component of B, B x . For ip > 0.5 the profiles are zero. 

negligible except in extreme cases. The solutions in this paper are therefore determined by 
the /photo boundary condition, a simple bipole, and the form of the governing equation. The 
subphotospheric line current modeled by J p hoto provides a convenient means of producing 
the bipolar boundary conditions. Alternative bipolar boundary conditions would produce 
similar results. 

The four plots in Figure [3] are labeled a, b, c and d to match the four curves in Figure [2j 
Curve c has a shallower decrease and curve d a steeper decrease in axial flux as one moves 
away from the magnetic axis than curve a, while curve b has the same decrease. Models 
a and b are very similar indeed. There are differences but they are too small to see in the 
plots. Model c is larger than a and b while model d is smaller. These effects are due to the 



Fig. 3. — Force-free equilibria of inverse-polarity prominence magnetic fields with f 2 given 
by the curve in Figure [2] with label a (top left), b (top right), c (bottom left) and d (bottom 
right). 

different quantities of axial magnetic flux confined close to the magnetic axes of the models. 
Figure H] shows the angle between the magnetic vector and the magnetic axis in each model. 
Here, the difference between models a and b is much more obvious. The magnetic angle of 
model b has a significantly sharper change at the edge of the flux rope than is the case for 
model a, and its angle is consistently smaller over much of the flux rope volume, indicating a 
more sheared field. In model d almost all of the angle changes are concentrated at the center 
and the edge of the flux rope. Here there is a concentration of highly sheared field at the 
center of the flux rope, outside which the angle has a plateau between 60° and 70°, finally 
increasing sharply to 90° at the edge. In the other models the distribution of angle changes 
is more even. The clumping of angle contours near the edge of each flux rope in models b, c 



Fig. 4. — Force-free equilibria of inverse-polarity prominence magnetic fields: angles between 
the field vectors and magnetic axes in the four models of Figure EJ The labels a, b, c and d 
match those in Figures [2] and [31 

and d is an effect of the total field strength there - small changes in f 2 cause large changes 
in the angle there. 

Chromospheric fibrils as observed in Ha are believed to delineate magnetic flux tube 
trajectories emanating from the photosphere. Over the neighborhood of an active region 
these fibrils follow an organized pattern: approaching the filament from either side, the 
fibrils become longer and aligned more parallel to the filament axis. Far from the filament 
the fibrils are nearly perpendicular to the filament, while at either side of the filament 
they are parallel. No fibrils close to the polarity inversion line are observed to arc over it, 
and the fibrils indicate that magnetic flux in the positive and negative magnetic regions is 
directed away from and towards the inversion line respectively. This is all contrary to what 
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is expected of a near-potential magnetic field and indicates the presence of significant electric 
currents in the vicinity of as well as within the filaments. Low & Hundhausen (1995) and 
Low & Zhang (2004) included photospheric magnetic shear in their models via their current- 
carrying fluxrope field. If this field is allowed to dip beneath the photosphere, the solution 
has a sheared photospheric signature similar to the sheared photospheric polarity inversion 
lines associated with prominences in observations. The same effect can be produced by our 
models by setting the photosphere to be just above z = 0.5. 

We also tried to find examples where the direction of the axial flux reverses. One im- 
mediate motivation for doing this is the observation by Lin et al. (1998) of a prominence 
in absorption against the disk in mid-latitude with axial field component of opposite signs 
on opposite sides of the photospheric neutral line below the prominence. Flux rope con- 
figurations resembling the reversed-field pinch (Freidberg 1987, Biskamp 1993) familiar to 
plasma laboratory physicists would exhibit precisely this pattern. In Taylor's (1986) theory 
of plasma relaxation to minimum energy states under conservation of magnetic helicity, there 
exists a helicity threshold above which the relaxed field must exhibit a reversal of its axial 
field across the tube cross-section, a result repeatedly confirmed by experiment. In the more 
complex, open-state environment of the solar atmosphere it is not clear if corresponding 
phenomena are to be expected. We were not able to obtain to obtain cases with significant 
reverse axial flux confined by a potential ambient field. The necessary decrease of f 2 (il>) 
near the center of the flux rope causes the contour surfaces to crowd near the magnetic axis 
and the consequent strength of the non-axial flux makes deviations of the magnetic vector 
from the y-z plane of more than 10° very difficult to obtain. One simple way to construct a 
model with reversing axial flux would be to include axial flux in the overlying bipolar arcade, 
of opposite sign to the axial flux in the flux rope. A sheared flux rope emerging beneath 
an arcade of opposite shear may be the likeliest scenario for the creation of structures with 
reversing axial flux in the open solar corona. We do not pursue this option here. 



4.2. Case in which temperature T = T(ip) is a flux function 

We now turn our attention to non-force-free MHS solutions. For the present we set 
f 2 (ip) to be as in Equation (TTOT) with n = 2. In the isothermal case, the hydrostatic scale 
height A = p/(pg) is constant and we may set t(z) = exp(— z/A). Now the condition 
p/p = constant gives pi(ip) = and equation ([8]) integrates to give 



P = Po(^)exp(-z/A), 



(11) 
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consistent with equation (j3J). 

In the more general case with T = T(tp), Equation (j7|) integrates to give 

p(ip, z) = p Q (ip) exp ( —r^-r ) , 




(12) 



where A(^) = ksT(ip) / '(fig). The Grad-Shafranov Equation ([6]) becomes 
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The isothermal case is T(ip) = constant or A(-0) = constant, which coincides with the 
case of Low (1975). The additional term for non-isothermal cases models the pressure gradi- 
ent introduced by temperature changes across field lines. In a system including very different 
temperatures, such as a cool prominence embedded in a much hotter coronal medium, this 
term may be important. 

The forms for po(ip) considered in this work are plotted in Figure [51 Curve a is the 
form used by Low & Hundhausen (1995) and Zhang & Low (2004): po(ip) = 3/10 — ip. For 
ip > 3/10, po = and the magnetic field is force- free there. Curve b is the graph of 



Curve b is a smooth alternative to curve a, sharing its value and derivative at ip = but 
having a double instead of single root at ip — 3/10. We adopt this form of po(ip) for the 
examples plotted in Figure El 

The examples shown in the top two panels of Figure [6] are solutions to the same equation 
with the same boundary conditions but with different values for the temperature. These first 
two examples have no axial magnetic flux and have smaller plasma (3 = 87rp/|B| 2 than the 
sheared examples plotted below them. The examples of the bottom pictures the plasma (3 
has maximum of about 1/10 or 1/5. The (3 parameter is infinite at the unsheared solutions' 
neutral points at the center of the flux ropes, but a typical value for (3 in these solutions is 
about 10. The Bernoulli integral along field lines is 



PoW = t^(3-10^) 2 (3 + 10^). 



(14) 




(15) 
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Fig. 5. — Examples of po(ip) functions applied in this paper. For ip > 0.3 the profiles are 
zero. 

where po is the density at z = 0. This expresses the fact that the thermal energy and grav- 
itational potential energy must add to give a conserved quantity the total energy. Clearly 
the thermal energy is directly proportional to the hydrostatic scale height and therefore the 
temperature. The top pictures of Figure M show that an increase of thermal energy of the 
plasma causes the thermal pressure to be more effective in pushing out against the confining 
forces of the magnetic field. The more even spread of mass across the concave-upward floor 
of the magnetic configuration imposes less stress on the magnetic field at y — and allows 
the prominence plasma to sit higher in the atmosphere than in cases with lower temperature. 
The effect of this on the characteristic scale length of both the plasma distribution and the 
magnetic structure itself are clearly visible in the top two pictures. 

The bottom two pictures show examples with colder plasma but with significant axial 
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magnetic flux density - they have axial flux profiles derived from that of model b of Sec- 
tion 14.11 while the magnitudes of their p profiles are similar to those of the top pictures. 
These stronger magnetic structures are able to accommodate cold plasma without significant 
changes in their magnetic fields. In reality prominences have typical hydrostatic scale height 
of a few hundred km while flux ropes are of order 10 Mm across. The bottom left picture 
shows the model closest to these parameter values. 

To demonstrate the flexibility of the solution method, we show an example with two 
separate prominence concentrations in a single flux rope. Curve c in Figure is, 

200000 / 3 V / 3 ,\ 2 /3 



»<*> = — [jo-") U + (16) 

This polynomial has double roots at ip = 3/10 and 3/20 and has the same value and derivative 
at ij} — as 3/10 — ?/>. The corresponding prominence model is plotted in Figured A narrow 
evacuated channel, corresponding to the first zero of Po(ip), separates a rounded blob of 
plasma above and a curved sheet below. 

Further complexity is introduced by an ambient atmosphere populated by a hot plasma 
in a non- force- free magnetic field. Figure M shows a model which exploits the non-isothermal 
aspect of Equation f|T3|) . A ^-dependent temperature profile is imposed so that a hot corona 
and a cool, dense prominence can be modeled together in a self-consistent way. The variation 
of the hydrostatic scale height A(ip), proportional to the temperature, with ip within the flux 
contour ip = 3/2 is described by 



AM = J; + ^ (17) 

This profile is graphed in the left panel of Figure [H The scale height therefore ranges 
from 3/20 to 3 between the magnetic axis and ip = 3/2. Outside ip = 3/2 the temperature 
has a constant value of 3. On the other hand, the plasma pressure is controlled by choosing 
Po(ip), which in this model is based on curve d of Figure El Curve d is similar to curve b but 
has derivative -3 at ij) = 0. The equation for curve d is 



PoW = ^(3-10^) 3 , (18) 

so that this curve has a triple root at ip = 3/10. This po(ip) is zero outside the contour 
ip = 0.3. In order to generate an ambient atmosphere we add a small number, 0.0005, to 
Po(ip) in all space. Therefore, outside a certain contour, ip = 1.5, the pressure contours are 
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planes parallel to the photosphere. Within the contour ip = 3/10 is a cool, dense prominence 
structure gravitationally stratified with scale height about 1/20 of that in the ambient corona. 
One result of this difference in stratification between the interior and the exterior of the 
flux rope is that the prominence plasma resides in a cavity: the portion of the flux rope 
not occupied by dense prominence plasma is more sparsely populated even than the corona. 
This is an essential part of the three-part structure often observed in prominences and CMEs 
(Low & Hundhausen 1995, Gibson et al. 2006). Low & Hundhausen (1995) emphasized and 
gave physical meaning to the three-part structure of the prominence-flux rope-helmet system 
often observed. For simplicity, the prominence is usually studied without accounting for the 
low-density cavity in which it is observed to be embedded. While the large-scale corona 
varies in structural complexity with the solar cycle, the basic three-part helmet structure 
seems to persist throughout this cycle. Here it is seen as a natural feature of a low-density 
ambient atmosphere containing a flux rope with a temperature minimum and a plasma 
pressure maximum at the flux rope center. 

4.3. Case in which entropy s = pj p 1 = s(ijj) is a flux function 

In the polytropic case, Equation ([7]) takes the form 



7 dp 
7 — 1 dz 



?/>=const 



(19) 



where p(ip, z) = s{jp)p(ip, z) 1 . In the special case of Low & Zhang (2004), po(if)) = po(ip)g/n 
for n = 7/(1 — 7). 



The Bernoulli integral along field lines is 



dp 7 

+ gz 



p 7-1 



(v__ 


p 




\P 


p 


z=o) 



+ gz = constant, (20) 



expressing the fact that the thermal energy and gravitational potential energy must add to 
give a conserved quantity. For 7 < 1, p/p — (p/p)\ z=Q must be > 0, i.e. the temperature 
increases with height while for 7 > 1, p/p — (p/p)\ z=0 must be < 0, i.e. the temperature 
decreases with height. Note that along field lines the temperature is a linear function of 
height, which must therefore cross zero and become negative at some height. Such solutions 
may still be used to model solar prominences since they are of finite vertical extent. In the 
separable case, the plasma pressure and temperature are given by 
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p{ip,z) 



7 - 1 (z - zb) 



7 

V (7 - 1) 



k, 



7 



A 



1 

7-1 



(21) 
(22) 



where A = p (ip)/p (ip)g is a constant length scale. The Grad-Shafranov Equation ([2D is 



dip 



+ 47T 



7 - 1 (2 - z ) 



7 



A 



dpoW 
dip 



(23) 



Note that the temperature is independent of ip in this separable case, since p and p have 
common -^-dependence, and the position of this zero-point is z = zq on all field lines. In order 
for the plasma pressure and temperature to be well-defined and positive in the domain of 
interest, the factor —(7 — l)(z — z ) must be positive throughout that domain. The gradient 
of the linear variation of T with altitude is determined entirely by the polytropic index 7. T 
is an increasing function of height for 7 < 1 and a decreasing function of height for 7 > 1. 
Thus, for example, Low & Zhang chose their z to be negative to keep this critical point out 
of the way beneath the base of the corona because they were working with values of 7 < 1. 
For cases with 7 < 1, Zq must be placed above the domain of interest and for 7 > 1 z 
must be placed beneath. If 7 > 1 this introduces difficulties for the construction of solutions 
in which all of space is populated by plasma, such as the one in Figure El but for isolated 
prominence plasma enhancements in empty atmospheres these solutions are appropriate. 



For the general case p 
temperature are given by 



Po(ip), the solution is not separable and the pressure and 



p(ip,z) 
T 



7 - 1\ (z- z ) 
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where 
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is the temperature as a function of tp at z — zo, and A(ip) = Po{i^)/ Po{ip)g is the hydrostatic 
scale height at z = Zq on each flux surface labeled ip. These are both zero in the separable 
case discussed above. The non-separable Grad-Shafranov Equation is 



In this general case T still varies linearly with z on each individual field line with gradient 
unaffected by T (i/j), but that T can now vary across field lines via T (ip). The effect of T (ijj) 
is therefore simply to shift the critical point z = z' = z + [7/(7 — l)]A(ip) vertically by a 
distance [7/(7 — l)]A(i/)) on each field line. This quantity can vary freely from field line to 
field line but its physical effect is indistinguishable from the constant critical point of the 
separable case z = z : fixing the temperature to be a given temperature at a given height 
fixes z' exactly as z is fixed in the separable case. 

Some illustrative polytropic hydrostatic ID solutions are shown in Figure [9) The location 
of a point z at a particular chosen temperature (in this case 7000 K, a typical temperature for 
prominence plasma) relative to the critical point Zq is graphed as a function of 7 in the top 
left picture. As anticipated earlier, for 7 < / > 1, zq is below/above the domain of interest. 
Since the temperature gradient is larger for values of 7 further from 1, the critical point must 
be positioned closer to the domain of interest for these values of 7 and far from the domain 
of interest for 7 close to 1. In the limiting case 7 = 1, the critical point asymptotes, crossing 
from —00 to 00. 

In Figure [9] the plasma parameters p, p and T for different values of 7 are forced to be 
equal at z — 0. Moreover, from the balance of pressure gradient and weight, the vertical 
derivatives of the pressures are also equal at z = so that all pressure curves are tangential 
there. Among the pressure curves is the 7=1 isothermal curve which has constant pressure 
scale height as discussed in the last subsection. Solutions with 7 > 1 are cooler than the 
7 = 1 solution in z > 0, have smaller scale heights and therefore have a steeper appearance 
there. On the other hand, solutions with 7 < 1 are hotter than the 7 = 1 solution in z > 0, 
have larger scale heights and therefore have a flatter appearance. Meanwhile, in z < the 
opposite is true: solutions with 7 > / < 1 are hotter/cooler and have flatter/steeper pressure 
curves. 

The density gradients, on the other hand, are not equal at z = since they are affected 
by the different temperature gradients of the different solutions. Density gradients are larger 
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for small 7 values than for larger 7 values since the plasma is cooler. They are also larger 
within a scale height (~ 400 km) of z = in the positive direction, where the temperature 
gradient causes the density to fall off more quickly for 7 < 1 and more slowly for 7 > 1 than 
in the 7=1 case. Beyond a scale height, the higher temperature of the 7 < 1 plasma causes 
these density curves to become flatter than the 7 = 1 curve and the lower temperature of 
the 7 > 1 plasma causes the density curves to become steeper than the 7 = 1 curve. 

Figure [TU] shows magnetohydrostatic 2D examples with a range of values of 7. The ID 
hydrostatic solutions of Figure and the 2D magnetohydrostatic solutions here are related in 
that the 2D solutions restricted to flux contours ip = constant are ID hydrostatic solutions 
of the type described above. The 2D solutions were all derived with the same boundary 
conditions as those presented in earlier sections and the differences between the plots in 
Figure (TU] are due only to the different values of 7 and the locations of z . In the 7 = 0.8 
case Zo must be below the flux rope while in the cases with 7 > 1 it must be above. In 
the four examples shown, zq is placed about a unit from the center of the flux rope. This 
means that for 7 > 1, the temperature increases as 7 increases, as can be seen in the plots. 
In these cases the hotter plasma is at the bottom of the flux rope while the hotter plasma 
is at the top in the 7 < 1 case. This makes the 7 > 1 plasma pressure distributions more 
bottom-heavy than the 7 < 1 case. 



4.4. A prominence with normal topology 

In the non-force-free models of this paper calculated in the Cartesian half-space z > 0, 
there are two current systems: one centered at the center of the flux rope associated with 
the axial magnetic flux and one associated with the body of plasma centered in the bottom 
half of the flux rope. From a distance, these current systems resemble line currents to lowest 
order, and their interaction determines the global field topology and far field structure of a 
solution. Thus the topology can be controlled via the signs and ratio of these currents. In 
inverse-topology cases the shear- and plasma-induced currents are flowing in the same axial 
direction. Normal-topology cases have additional complexity because the axially-directed 
currents associated with the shear and the plasma forces are flowing in opposite directions. 

An example is shown in Figure [TT1 FINESSE is less suited to the calculation of normal 
prominences than of inverse prominences because it is designed to find equilibria with ip = 
at the magnetic axis and ip > everywhere else, with maximum on the boundary. The flux 
function of a normal prominence has a local minimum at the magnetic axis but, since the 
overlying bipolar arcade and the lower half of the flux rope must have flux traveling in the 
same horizontal direction, the flux function must have a saddle point above the flux rope, 
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above which the flux function must decrease with height (see Figured]). The example in 
Figure [11] has such a saddle point far above the flux rope outside the field of view of the plot. 
It is not yet clear how to incorporate this and other topological complications in a controlled 
way. 

5. Discussion 

A major obstacle in the way of understanding solar activity is the difficulty of capturing 
in models the complexity of the plasma dynamics and of the magnetic field structure. For 
this reason, most recent efforts to model the solar atmosphere have attempted to recreate 
as far as practically possible the full 3D geometrical complexity of the physical parameters 
(e.g. Amari et al. 2003a, 2003b, Roussev et al. 2003, Wiegelmann & Neukirch 2006). Here 
we have adopted the alternative approach, following Low & Hundhausen (1995) and Low 
& Zhang (2004), of studying a simple generic physical system: a single solar prominence 
plasma enhancement suspended in a near-potential coronal magnetic field. 

We have given new numerical magnetohydrostatic solutions describing the gravitation- 
ally stratified, bulk equilibrium of cool, dense prominence plasma embedded in the near- 
potential coronal field. These solutions are calculated using the FINESSE magnetohydro- 
dynamics equilibrium solver and describe the morphologies of magnetic field distributions 
in and around prominences and the cool prominence plasma that these fields support. The 
new solutions were not accessible by previous analytical techniques. The numerical method 
allows us to prescribe in a controlled way the temperature or entropy as a function of the 
magnetic flux function, enabling flexibility of choice of the physical parameter distributions. 
We focussed on new solutions with a range of values of the temperature, the magnetic shear, 
the polytropic index and with large temperature variations perpendicular to the magnetic 
field. The axial component of the magnetic field gave demonstrably increased structural 
integrity to that field. With temperature a flux function and a low-density ambient atmo- 
sphere, a temperature minimum and a plasma pressure maximum at the flux rope center were 
found to produce a relatively evacuated cavity within which the cool plasma was embedded. 
Such a cavity has been seen many times in observations. For solutions with entropy a flux 
function, the polytropic index and the temperature distribution were found to be related in 
a simple way for separable cases. 

The stability properties of these new equilibria can be determined in the form of a 
full resistive MHD spectrum by solving the linearized MHD equations using a companion 
hyperbolic stability solver, PHOENIX (Blokland et al. 2007a). Because of the centrality of 
prominences and their fields to the most geoeffective space weather phenomenon, the CME, 



-23- 



precise knowledge of these stability properties is very desirable. Measurement of the spectrum 
of MHD waves, called MHD spectroscopy in analogy with quantum mechanical spectroscopy, 
is still in its infancy but may well lead to a firm knowledge of the internal characteristics 
of plasmas structures (Goedbloed 2002). MHD spectroscopy entails a separate study of 
the nonlinear static equilibrium configuration on the one hand and the various linear wave 
structures that can occur on the other. The first study is the subject of this paper and the 
spectroscopy will be treated in a sequel. The new equilibria also serve as suitable starting 
points for time-dependent MHD simulations using, e.g., the Versatile Advection Code (VAC). 

A 3D extension of this work is possible via the study of helically symmetric MHS 
configurations, which are also governed by a Grad-Shafranov equation with a linear elliptic 
operator. Bearing in mind the difficulty of calculating 3D MHS equilibria of mixed elliptic- 
hyperbolic type it is tempting to settle for a simple option: that of abandoning the direct 
calculation of equilibria and their linear perturbation spectra and exclusively working on 
time-dependent MHD simulations instead. This would be an inferior option since it passes 
up the precision and detail offered by the equilibrium and spectral theory. Optimally one 
would hope to exploit both approaches in tandem. 

We thank the referee for helpful and constructive comments. This work was conducted 
while GP was a participant in the National Aeronautics and Space Administration (NASA) 
Postdoctoral Program at Goddard Space Flight Center, and was based at National Solar 
Observatory, Tucson. GP thanks RK for kind hospitality during a research visit to Rijn- 
huizen. 



Suppose we seek only separable solutions to Equation (jSJ). In this case solutions of 
Equation (CO) for p and p must have identical ^-dependence. Let 



where t p (z) describe the field-aligned variation of the pressure with height. Then Equation 
takes the form 



A. Appendix 



A.l. Numerical solution of the problem 



(Al) 



V 2 ij> = -A[Fty) - Br p (z)Il(il>)], 
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where 
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and A and B are constants: B is related to the plasma (3 = 87rp/|B| and is specified at input, 
while A is to be determined by the numerical algorithm as an eigenvalue of the problem. 



The weak form of equation flA2[) is 



or 



f x V 2 i/jdV = -A [ X [FW ~ Br p (z)U^)]dV 
Jv Jv 

Vx ■ Vtj)dV = A [ x[F{ip) - BT p {z)H{ip))dV + [ xdn?pdS. 

JV JdV 



(A5) 



(A6) 



for arbitrary test functions X- Since this is true for any x, we may transform x $X m the 
equation above to obtain the variational form. 

We represent y, z, ip, x by bicubic isoparametric elements defined in a rectangular 
coordinate system (s,t) G [— 1, l] 2 . Within a cell, any function f(y,z) has finite element 
representation 



f(y(s,t),z(s,t)) = 



df df 
H o(s, t)f(y , zq) + H 10 (s, t) — (y Q , z ) + H 01 (s, ^^U/o, z o) 



d 2 f 



+H n (s,t)——(y ,z ) 
oyoz 



(A7) 



where H 00 , H 10 , . . . are the usual Hermite bicubic polynomials. The summation runs over 
corner points (sq, to) of the cell and yo = x(so,t ), zq = y(so,t ). Also, ip, y and z are 
approximated by the same locally defined interpolating functions 



dtp dtp d 2 ip 

Hoa(s,t)ip(s ,to) + H w (s,t)—(s ,t ) + H 01 (s,t)—(s ,t Q ) + H n (s,t)-^-^-(so,t Q ) 

(A8) 
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problem of the form K^x^J 1 ^ = b\ . Here 



etc. The Picard iteration of the variational form of equation (1A6|) leads to an iterative matrix 

L ij x j 



K lj= ! f VHi(s,t) ■VH j (s,t)Jdsdt, (A9) 
Jo Jo 

the 2( n+1 ) are the coefficients of the interpolating functions of the finite elements and 

b\ n) = a[ [ H^F^) - BT p {z)U{^)]Jdsdt, (A10) 
Jo Jo 



where J = |^fj and Hi are the finite elements, given single subscripts for simplicity. 



A. 2. Scaling 

We now introduce scale parameters which may be used to normalize the physical quan- 
tities. If Bo is a typical magnetic field strength, all quantities can be made dimensionless as 
follows: 



y = y/H, (Ail) 

z = z/H, (A12) 

/ = Sf/B , (A13) 

p = 5 2 Att P /B^ (A14) 

Here 5 = is a dimensionless parameter that enters when we choose to normalize the 

magnetic flux ip = ip/ipo where ipo is a typical ip value. In practice we only work with 
dimensionless quantities and drop the bars without confusion. 
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Fig. 6. — Magnetohydrostatic isothermal equilibria of inverse-polarity prominence magnetic 
fields: unsheared examples with A = 0.3 (top left) and 1.0 (top right), and sheared examples 
with A = 0.03 (bottom left) and 0.1 (bottom right). The flux contours, which coincide with 
magnetic flux trajectories projected onto the y-z plane, are represented by solid lines. The 
plasma pressure is represented by the shading, with dark indicating high pressure. 
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Fig. 7. — The flux contours and plasma pressure of a magnetohydrostatic isothermal equilib- 
rium with an inverse-polarity prominence magnetic flux rope containing a double prominence 
enhancement. 




Fig. 8. — A cool, dense prominence plasma enhancement suspended in a hot, sparse corona. 
The left picture shows a highly non- isothermal scale height profile. The right picture shows 
the flux contours and the plasma pressure. Note the evacuated cavity within which the dense 
plasma is embedded. 
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Fig. 9. — Hydrostatic ID polytropic equilibrium: for a fixed temperature at z — 0, the 
position of the critical point z = z is graphed as a function of the polytropic index 7 (top 
left). The resulting temperature (top right), pressure (bottom left) and density (bottom 
right) are shown as functions of height. 
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Fig. 10. — Magnetohydrostatic polytropic equilibrium: 7 = 0.8 (top left), 1.1 (top right), 
1.4 (bottom left) and 5/3 (bottom right). The critical point z is located at z = for the 
7 = 0.8 case and at z = 2 for the 7 > 1 cases. 
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Fig. 11. — Magnetohydrostatic isothermal equilibrium with normal topology. 



